A review of Monte Carlo simulations of polymers with PERM 
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In tliis review, we describe applications of the pruned-enriclied Rosenbluth method (PERM), a 
sequential Monte Carlo algorithm with resampling, to various problems in polymer physics. PERM 
produces samples according to any given prescribed weight distribution, by growing configurations 
step by step with controlled bias, and correcting "bad" configurations by "population control". The 
latter is implemented, in contrast to other population based algorithms like e.g. genetic algorithms, 
by depth-first recursion which avoids storing all members of the population at the same time in 
computer memory. The problems we discuss all concern single polymers (with one exception), but 
under various conditions: Homopolymers in good solvents and at the Q point, semi-stiff polymers, 
polymers in confining geometries, stretched polymers undergoing a forced globule-linear transition, 
star polymers, bottle brushes, lattice animals as a model for randomly branched polymers, DNA 
melting, and finally - as the only system at low temperatures, lattice heteropolymers as simple 
models for protein folding. PERM is for some of these problems the method of choice, but it can 
also fail. We discuss how to recognize when a result is reliable, and we discuss also some types of 
bias that can be crucial in guiding the growth into the right directions. 
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I. INTRODUCTION 

Research in the field of polymer physics has grown vig- 
orously since the 1950s [iHj]- Recent developments in 
the techniques for the tools of atomic force microscopy 
(AFM) [sll, in fabrication of nanoscale devices and in 
single-chain manipulation techniques |6|-l8| open possibil- 
ities for a broad range of applications in physical chem- 
istry, biotechnology and material science. During this 
time, much effort has also been put into studying the 
statistical properties of polymers by computer simula- 
tions [3, [lO]- Indeed, due to the richness of the ob- 
served phenomena and the non-triviality of the problems 
involved, polymer physics has from the very beginning 
served as a playgr ound for developing novel Monte Carlo 
strategies |lll - [l3 |. These strategies depend strongly on 
the problems one is interested in: Linear versus branched 
polymers, dilute versus dense systems, scaling laws ver- 
sus detailed material properties, classical versus quantum 
mechanical problems, implicit versus explicit treatment 
of solvent, etc. 

In this review we shall only deal with one class 
of algorithms, the pruned-enriched Rosenbluth method 
(PERM) d4]. So far it has been used for classical 
physics only, although closely related methods have also 
been used since long ago for quantum mechanical simu- 
lations [ia| . Although it is not a panacea and fails miser- 
ably in many problems, it still found applications to sev- 
eral of the above dichotoma, and in some cases it beats 
the (presently known) competitors by huge margins. 

In the following we shall mostly be concerned with sin- 
gle unbranched molecules moving freely in a dilute sol- 
vent. Later we will also consider branched polymers and 
polymers attached to surfaces. The basic characteristics 
of linear polymer chains depend on the solvent condi- 
tions. At high temperatures or in good solvents repulsive 
interactions (the excluded volume effect) and entropic ef- 



fects dominate the conformation, and the polymer chain 
tends to swell to a random coil. At low temperatures or 
in poor solvents, however, attractive interactions between 
monomers dominate the conformation and the polymer 
chain tends to collapse and form a compact dense globule. 
The coil-globule transition point is called the 0-point. 
Based on field theory Q , the behavior of polymer chains 
in good solvents is well understood. In the thermody- 
namic limit (as the chain length N — >■ cx)), the partition 
function scales as 



fi^^'N^-^ at T > Te 



(1) 



where fioo is the critical fugacity and 7 is the entropic 
exponent related to the topology. Below the 0-point, a 
collapsed polymer can essentially be viewed as a liquid 
droplet. According to the Lifshitz mean field theory [2, 
IJ] , a surface term should be included in the partition sum 
as 



Z^a^b^^N^- 



at T < Te 



(2) 



with s = {d — l)/d and b > 1. 

Generally speaking, the thermodynamic limit of a 
polymer system coincides with the limit when the chain 
length N tends to infinity. For conventional Monte Carlo 
(MC) methods such as the Metropolis algorithm, one can 
only simulate moderately large systems, the maximal fea- 
sible values of N depending on the temperature and on 
the degree of reality of the model. Going from simple 
lattice-based models at high temperatures to models with 
realistic interactions and further to folded proteins with 
explicitly included solvent, iVmax niight decrease from 10^ 
to ^ 100. If one is interested mostly in scaling laws (as 
we shall be), one simulates at several values of N and 
uses finite-size scaling (FSS) to extrapolate the behavior 
of the considered thermodynamic quantities to iV — > 00. 
Rather often, either very large finite-size effects have to 



be considered or it is too difficult to reach equilibrium 
states or to produce sufficiently many independent con- 
figurations. For some problenis (not for all!), it was a big 
breakthrough when (PERM) [IJ) 11S43 ^^^ proposed in 
1997. It is particularly efficient for temperatures near the 
Q collapse, where chains of length up to A^ = 1, 000, 000 
could be sampled with high statistics, and it was con- 
firmed unambiguously that the Q collapse is a tricritical 
phenomenon with upper critical dimension dc = 3 [3|. 
Since then, many other applications have also been made. 
Many other applications have also been made successfully 
by PERM [ig , which provide in some cases a deep under- 
standing on the scaling behavior of polymer chains under 
different solvent conditions, geometrical confinements, on 
the phase transition behavior of polymer chains adsorbed 
onto a wall, on polymers stretched by a force, etc. 

In the next section we give a detailed description of 
the basic algorithm. This algorithm can be made sub- 
stantially more efficient by a suitable bias in the growth 
direction, and two biases (including 'Markovian anticipa- 
tion') are discussed in Sect. 3. Applications are treated 
in Sects. 4 (0-polymers), 5 (stretched polymers in poor 
solvents), 6 (semifiexible polymers), 7 (polymers in con- 
fining geometries), 8 (branched polymers with fixed tree 
topologies), 9 (lattice animals), 10 (protein folding), and 
11 (DNA melting). Finally the paper concludes with a 
summary in Sect. 12. 



II. ALGORITHM: PRUNED-ENRICHED 
ROSENBLUTH METHOD (PERM) 

In statistical thermodynamics, the partition function 
for a canonical ensemble in thermal equilibrium is defined 
by 



ZW) - Y. Q(") = E exp(-/3i?(a)) 



(3) 



here /3 — l/ksT, E{a) is the corresponding energy for 
the a^^ configuration, Q{a)/Z is the Gibbs-Boltzmann 
distribution, and Q{a) = exp(— /3i^^(a)) is normally 
called the Boltzmann weight. If each configuration is 
repeatedly and independently chosen according to a ran- 
domly chosen probability p{a) (a bias), the partition sum 
is rewritten as 



Z = lini Z 



(4) 



where M is the number of trials and 

^ = T7 E Q(«)/P(") - TF E ^(") ' (5) 



Q— 1 a.=X 

with modified weights 

W{a) = Q[a)/p{a). 



(6) 



If we use p[a) ex exp(— /3£^(a)) {Gibbs sampling}, each 
contribution to Zm has the same weight, which is an 



example of the so called 'importance sampling'. The es- 
timate for any observable A is given by 

(A) . hm (A)m = hm i:^ii/(-)^(") . (7) 

In general, we expect that statistical fluctuations of {A) m 
are minimal, at given M , if we use importance sampling 
and if all trials are independent. In general this is infeasi- 
ble. The Metropolis method achieves perfect importance 
sampling at the cost of highly correlated trials. PERM 
tries, with a completely different strategy, at a compro- 
mise between importance sampling and independence. 

Things are best illustrated by a linear chain of A^ -f 1 
monomers in an implicit solvent, modeled by an inter- 
acting self-avoiding walk (ISAW) of N steps on a sim- 
ple (hyper-)cubic lattice of dimension d. The interac- 
tions in this model are (i) the chain connectivity which 
enforces that adjacent monomers sit on adjacent lattice 
sites; (ii) self-avoidance that excludes configurations in 
which the same lattice site is occupied by two or more 
monomers; and attractive interactions (energies e < 0) 
between non-bonded monomers occupying neighboring 
lattice sites. Writing q = exp(— /?e) for the Boltzmann 
factor, the partition sum is 



ZN{q) 



walks 



(8) 



where m denotes the total number of non-bonded nearest 
neighbor pairs. The solvent quality is varied by changing 
the temperature T. 

In the original Rosenbluth-Roscnbluth (RR) 
method [llj, polymer chains are built like random 
walks by adding one monomer at each step. At the 
0*'^ step, the first monomer is placed at an arbitrary 
lattice site. For this "chain" of length A'' = 0, the 
weight is trivially Wq = 1. For the first step one has 
2d possibilities and no interactions yet, giving Wi — 2d. 
For subsequent steps one has to take self-avoidance into 
account. When a monomer is added to a chain of length 
A^ — 1, one scans the neighborhood of the chain end 
to identify the free sites on which a monomer could be 
added. If there are rifiee > 1 free neighbors, the next 
step is chosen uniformly among them, while the walk 
is killed if rifree = ("attrition"). After this step the 
weight Wn is updated by multiplying Wn-i by 



w„ 



q'^-fPn 



(9) 



where p„ = l/rifroe and r7i„ is the number of neighbors of 
the new site already occupied by non-bonded monomers 
(notice that m = X)n=o"^")- Therefore, after A^ steps 
the total weight is 



Wn = wnWn- 



WnWN-1 



.Wo 



N 



(10) 
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FIG. 1. Schematic drawings of building a SAW from the 0**^ step to the 4"^ step and the associated weight at each step. At the 
O' step, we set the weight Wo — wo — 1, and the probability po = 1. The Rosenbluth bias is used here such that pn = l/wfroe 
at each step, so the total weight W„ — Tl^^Zonhcc- 



When the chain length A^ becomes very large, the method 
fails for two reasons: First of all, attrition implies that 
only an exponentially small fraction of trials survive and 
give any contribution at all. Secondly, as the weight fac- 
tors Wn are weakly correlated random variables, the full 
weight Wn will show huge fluctuations. Thus the surviv- 
ing configurations will finally be dominated by a single 
configuration, demonstrating a dramatic lack of impor- 
tance sampling. 

PERM [IJj was invented to overcome this shortcom- 
ing of the RR method. The main spirit of PERM is as 
follows, 

• Polymer chains are built like random walks by 
adding one monomer at each step. 

• A Rosenbluth-like bias is used for choosing one of 
the nearest neighbor free sites for the next step of 
the walk, but a wide range of probability distribu- 
tions can be used depending on the specific prob- 
lem at hand, which will be discussed in detail in 
the following sections. 

• In order to overcome attrition and to reduce the 
fluctuations of W„, one uses "population control". 
This is achieved by pruning some low-weight con- 
figurations and cloning ("enriching" [IJ]) all those 
with high weight, as the chain grows. To define 
'low' and 'high' weights, one uses two thresholds 
W+ and W~ . If at any step n the current weight 
Wn according to p^ would be > W^ , we make 
k additional copies (typically A; = 1) of the cur- 
rent configuration and give each copy a weight 
Wn — WnWn-i/{k + 1). On thc other hand, if 
([To)) would give Wn < Wn , we call a random num- 
ber r G [0,1]. If r < 1/2 we kill the configura- 
tion. Otherwise, we keep it and double its weight, 
Wn = 2wnWn-i- It is casy to see that pruning and 
cloning leave all averages unchanged. It improves 
importance sampling enormously, but it also leads 
to correlated trials. 

For most problems the choice of the thresholds W^ 
and Wn is unproblematic, and they can be chosen 
simply as constant multiples of the the current es- 
timate of the partition sum given by ^, 



w: 



C+Zn 



and 



w- 



C-Z„ 



(11) 



were C+ and C_ are constants of order unity. A 
good choice for the ratio between C+ and C_ is 
found to be C+/C- ^ 10 in most cases. If (fTTj) 
does not lead to good results, chances are that the 
method would not work with any other choice ei- 
ther. If the method works well, PTj) gives sam- 
ples where the total number of length n configu- 
rations is independent of n, i.e. attrition is com- 
pletely eliminated and pruning & cloning compen- 
sate each other exactly (up to statistical fluctua- 
tions), for large n. For the first trials (when there 
is not yet any estimate Zn)^ we choose normally 
Wn = and W+ = oo (a very large number like 
10^"*^), which gives the original RR method. 

• The copies made in the enrichments are placed on 
a stack, and a depth-first implementation is used 
for the chain growth: At each time one deals with 
only a single configuration until a chain has either 
grown to the maximum length N or has been killed 
due to attrition. If the first happens or if the stack 
is empty, a new trial is started. Otherwise, the 
configuration on top of the stack is popped and 
the simulation continues. This is most easily im- 
plemented by recursive function calls. Since only 
a single configuration has to be remembered dur- 
ing the run, this requires much less memory than a 
breadth-first implementation that uses an explicit 
"population" of many configurations, as it is tradi- 
tionally used e.g. in genetic algorithms. 

• As we said, configurations obtained from different 
clones of the same ancestor will not be uncorre- 
lated. The set of all such configurations is called a 
"tour" . Different tours are uncorrelated. Depend- 
ing on the amount of cloning/pruning, however, the 
correlations within a tour could be so strong as to 
render the method obsolete. In that case the dis- 
tributions P(ln(yV')) of logarithms of tour weights 
W is very broad, so that we are basically back to 
the problem of the RR method (with single tri- 
als replaced by tours): Averages might be domi- 
nated, in extreme cases, by a single tour. For check- 
ing against this, we simply look at P(ln(yV)) (see 
Fig. [2]), and compare it with the weighted distri- 
bution WPilnW). If yVP(lnW) has its maximum 
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FIG. 2. Histograms of logarithms of tour weights P(ln W) 
normalized as tours per bin, and weighted histograms 
W^P(ln W) are shown as indicated. Weights W are only fixed 
up to a /3-dependent multiplicative constant. The simulation 
shown in panel (a) is reliable, while that in panel (b) is not. 
Adapted from Ref. [3. 



at a value of InW where the distribution P{hiW) 
is well sampled, we are on the safe side. If not, 
then the results can still be correct, but we have no 
guarantee for it. An illustration of these two cases 
is shown in Fig. [5] [l^. Fig ^a) shows that the 
sampling is sufficient and the statistical weight dis- 
tribution is reliable, but Fig. [2jb) shows the oppo- 
site situation where the result might be completely 
wrong. 



III. BIASED GROWTH 

An important aspect of the method is that in general, 
for high efficiency, one should choose judiciously a bias 
in the growth, in order to reduce as much as possible the 
fluctuations of the weight factors Wn ■ The optimal choice 
of bias is often a result of trial and error, as there exists 
no general theory for it. The two choices discussed in the 
following subsections are often useful, but by no means 
in all cases - and other choices may be useful in other 



applications. 

One aspect of PERM that often decides the success 
or failure is that any bias that improves the growth at 
an intermediate stage should also be helpful later, i.e. it 
should not lead the growth into a dead end. One ap- 
plication where this is violated dramatically is e.g. the 
problem of random walkers in a medium with randomly 
placed traps (the "Wiener sausage" problem, leading to 
the famous Donsker-Varadhan stretched exponential sur- 
vival probability [20]). In this problem walkers should, to 
maximize their survival chance at very long times, stay 
very close to their starting point. On the other hand, 
for short times the path integral (partition sum) is dom- 
inated by walkers who venture out to explore a larger 
area, even if that might mean they get killed by a trap. 
Since this system can be mapped onto a polymer prob- 
lem, one can apply PERM to it [IJ. These PERM simu- 
lations gave indeed the first unambiguous numerical ver- 
ification of the Donsker-Varadhan law, nevertheless they 
completely failed for very long times, because both bias 
and population control conspired to "mislead" the walk- 
ers dll to venture too far out. 



A. Global Directional Bias 

Assume you want to simulate a polymer whose one end 
is held fixed at x = 0, and the other end is pulled away 
by a constant force F. In Sect. |V] we shall discuss in de- 
tail the case of a poor solvent where the stretching might 
unfold the dense globule into which the unstretched poly- 
mer would collapse. Here we just discuss qualitatively a 
polymer in a good solvent, i.e. a stretched SAW. 

This system could of course be simulated by an un- 
biased SAW, and the stretching could be taken into ac- 
count by reweighing each obtained configurations with a 
Boltzmann weight oc exp(— /3x • F). But this would be 
extremely inefficient for large F, since weights would be 
very uneven, and "correct" configurations would be very 
rare and would have very high weight. 

A much better strategy is to use a bias in the direction 
of the next step of the walk in the direction of F. The 
amount of the optimal bias cannot be determined a priori, 
but depends also on the excluded volume effect which 
helps to push the end further away in the direction of 
the bias. We do not show any data here, but we just 
mention that the simulations get easier with increasing 
F, since the walk resembles more and more an ordinary 
biased walk in this limit, and pruning/cloning events get 
more and more rare. 



B. PERM with fc-step Markovian anticipation 

A less trivial bias is suggested by the fact that a grow- 
ing polymer will predominantly grow away from the al- 
ready existing part of the chain. This could be modeled 
crudely by determining the center of mass of that part. 



and biasing the growth away from it. A better strategy 
is to learn on the fly how a typical short existing chain 
(of k monomers, say) would bias the further growth in 
detail, and to remember at any time the previous k steps. 
This is called Markovian anticipation [IS . 1 2 2:^24] . 

The crucial point of the fc-step Markovian anticipation 
is that the (fc + 1)"^ step of walk is biased by the history 
of the previous k steps, i.e., the bias depends on the last 
k steps. Let's consider the general case of a walk on a 
d-dimensional hypercubic lattice. At each step ?, a walk 
can move towards to one of the 2d directions denoted by 
Si = 0, . . ., 2(i — 1. All possible configurations of {k + 1) 
steps {i — —k, —k + 1, . . ., —1, and 0), which are in total 
(2d)'^+^ configurations, are labelled by 



S = (s_ 



S-l, So) = (s. So) 



(12) 



here s and sq denote the configurations of the previous 
k steps and the (fc + 1)"^ step. Either during a separate 
auxiliary run or during the first part of a long run we 
build a histogram i/,„(S) with (2^)*^+^ entries. For any 
S, the value of Hm{S) is the sum of all contributions to 
Zn+m of configurations that had coincided with S during 
the steps n — k, n — k + 1, . . ., and n, summed over all 
n in some suitable range excluding transients. Typical 
values for 3-d SAWs might be fc = 10, to = 100, n > 300. 
Then _ff,„(S)/iJo(S) measures how successful configura- 
tions ending with S were in contributing to the partition 
sum m steps later. The bias in fc-step Markovian antici- 
pation for the next step is thus defined by 



P(so|s) 



H„-,{s,so)/Ho{s,sq) 



4 = 



Hm{s,s'Q)/Ho{s,s'Q) 



(13) 
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FIG. 3. The mean square end-to-end distance R% plotted 
against A^ in a log-log scale (a) and against 1/lnA^ in the 



normal scale (b) 



R^/N ex 1 - 37/(363 In TV) is indicated 



by the straight line. Adapted from Ref. [14 



IV. e-POLYMERS 

The first application of PERM was to 0-polymers in 
three dimensions jljj. As we said, the upper critical di- 
mension for the Q collapse is d = 3, whence we expect 
ordinary random walk behavior with logarithmic correc- 
tions. These corrections have been calculated to lead- 
ing 25] and next-to-leading [26| orders. The experimen- 
tal verification of these corrections is highly non-trivial, 
because one has to use extremely diluted solutions in or- 
der to avoid coagulation of different chains, and thus the 
signals are very weak. Nevertheless, they have been ob- 
served in small-angle neutron scattering 27| . 



A. A Single 0-Polymer 

It is for this problem that PERM shows the biggest 
improvement over all other Monte Carlo methods. The 
reason is that at the Q point entropic and energetic 
(Boltzmann-) effects cancel exactly in the limit N —>■ oo. 
For finite TV they do not cancel exactly (this gives rise 
to the logarithmic corrections), but it is still true that 



the weight factors w„ are very close to 1. Thus hardly 
any pruning/cloning is needed, and to a first approxima- 
tion the simulation reduces simply to a straightforward 
simulation of random walks with small weight correc- 
tions. Full PERM simulations for very long chains (the 
longest chains in [1J| had A^ = 10^) do require in average 
one pruning/cloning step for every 2000 ordinary random 
walk steps. Therefore, in chain length n the algorithm 
effectively performs a random walk with diffusion coeffi- 
cient D « 2000. Asymptotically for A — > oo the algo- 
rithm still needs 0{N'^) steps to create one independent 
configuration of full length, but the coefficient is tiny. 

Indeed, since a growing polymer with endpoint in a 
locally denser region might feel an elevated Boltzmann 
factor at step n, but feels the compensating entropic dis- 
advantage only one step later, the optimal algorithm that 
produced these results was a slight modification of the 
algorithm described in the previous section, where the 
population control was based on a modified weight with 
incremental weight factors 



w„ 



^IV 



'n+l- 



(14) 



instead of ([9]). Results are shown in Fig. [3] The- 
ory [25i] predicts leading logarithmic corrections to be 
Rli/N ex 1 - 37/(363 In A^), which would be a straight 
line in Fig. [3{b) with very small negative slope. Com- 
pared to that, the corrections to random walk behav- 
ior seen in Fig. El^b) are much larger, although they are 
clearly smaller than one would expect for, say, a power 
law correction. It was indeed shown in [26| that the next- 
to-leading term increases the deviation from mean field 
behavior and improves thus the agreement between the- 
ory and simulation, but a fully quantitative verification 
remains elusive. 

Far below the Tq, PERM becomes inefficient, and it 
is instructive to see why: In strongly collapsed globules, 
polymer configurations are locally similar to those in a 
dense melt, and are well approximated by simple random 
walks without any correlations Q . But this implies that a 
collapsed chain with N = 1000 has a configuration that is 
completely different from the first 1000 monomers of, say, 
a collapsed chain of 8000 monomers. The former would 
form a compact globule, while the latter would form a 
rather dilute structure. Thus, similar to the problem 
discussed at the end of the last section, bias and pop- 
ulation control during the early stages of growth would 
be completely misleading as far as late stages of growth 
are concerned. Otherwise said, by effectively disallowing 
configurations that are initially dilute and fill the inte- 
rior only during the later growth, the entropy is severely 
underestimated. 



B. Unmixing Transition of Semidilute Solutions of 
very long Polymers 

Let us now consider a semidilute solution of polymers 
of common length iV slightly below the Tq temperature. 
The "unmixing" transition at which these polymers coag- 
ulate and phase separate from the solute is, for any finite 
chain length iV, in the Ising universality class [2^|. As 
N —i' oo, the transition temperature Tc should approach 
Tq from below. Since the Ising model has upper criti- 
cal dimension dc = 4, but the 0-point has upper critical 
dc = 3, all critical exponents referring to collective prop- 
erties (correlation length, specific heat) should be that 
of the Ising model, while properties characterizing the 
iV-dependence (e.g. radii of gyration, critical concentra- 
tion, Tq — Tc) should be mean field like with logarithmic 
corrections. In particular, the monomer density at the 
critical point should scale as 



$n 



7V-l/2^ 



(15) 



up to logarithms of N. 

A long standing problem in the 1990's was that all 
experiments showed <l>c ~ N^'^" with Xc ~ 0.38 ± 0.01 
[28| , which was considered as incompatible with theory - 
in particular, since experimenters viewed any prediction 
of logarithmic corrections with great skepticism. 



PERM can be easily modified for multi-chain systems, 
simply by placing the first monomer of a new chain not 
near the end of the last chain, and by applying the correct 
combinatorial factors that take into account the identi- 
ties of different chains [29j. Such simulations are very 
inefficient for short chains, since then Tc <C Tq, but they 
become more and more efficient as A^ — >■ oo . They showed 
clearly that the deviations from p3|) are not due to a 
different critical exponent, as was believed at this time, 
but due to logarithmic corrections 2a]- These are much 
larger than predicted by theory [30[, but this was to be 
expected in view of the results for single isolated chains. 



V. STRETCHING COLLAPSED POLYMERS IN 
A POOR SOLVENT 

As a collapsed polymer chain of chain length N is 
stretched by an external force under poor solvent con- 
ditions, one observes from a collapsed globule phase to a 
stretched phase, as the stretching force is increased be- 
yond a critical value [3ll |. This phase transition is first 
order in d = 3 dimensions, as is also suggested by the 
analogy of the Rayleigh instability of a falling stream of 
fluid, but it seems to be second order in d = 2 [31]. Here 
we shall only discuss the 3-d case. 

This is modeled as a biased interacting self-avoiding 
random walk (BISAW) on a simple cubic lattice in three 
dimensions. Assuming that a chain is stretched in the 
x-direction by the stretching force F — Fcx {fix is the 
unit vector in the x-direction), an additional bias term 
6^ is incorporated into the partition sum given by ([8|). 
where h — exp{/3aF) is the stretching factor (a is the 
lattice constant) and x is the distance (in units of lattice 
constants) between the two end points of the chain in the 
direction of F. The partition sum is therefore 



ZNiq,b)= J2 O" 



(16) 



alks 



The poor solvent condition is indicated hy q > qq where 
1.3087(3) 14]. According to the scal- 



qe 



^ p-t/kTe 



ing law 121), in the thermodynamic limit iV — >■ oo, the 
partition sum for polymers in a poor solvent scales as 

- In ZN{q, 5 - 1) « fioo{q)N + d{q)N^/^ - (7 - 1) IniV 

(17) 
with /ioo being the chemical potential per monomer in an 
infinite chain, and a is related to the surface tension a. 
Choosing q—l.b which is deep in the collapsed region, 
we performed simulations of BISAW with PERM. In or- 
der to improve the efficiency, each step of a walk is guided 
to the stretching direction with a higher probability. The 
n'^ step of walk (adding the (n-|-l)"^ monomer) is toward 
one of the free nearest neighbor sites of the v}^ monomer 
in the parallel, antiparallel, and transverse direction to 
F with probability: p^ : p^ : p± = \jh : \/l/b : 1. Thus 



we have 



(a) 



Pi = 



Z^allo 



if the step of the walk toward to 
the i — direction is forbidden 

otherwise 



(18) 
The corresponding weight factor at the n*'* step is then 



lAi 



Pi 



(19) 



where rra„ is the number of non-bonded nearest neighbor- 
ing pairs of the (n + 1)*^ 
the displacement (r„+i - 

n**^ monomers is in the direction perpendicular, parallel, 
and antiparallel to F, respectively. The total weight of a 
chain of length n is then 



monomer. Axi = 0, 1 or —1 if 
r„) between the (n + 1)*'^ and 



Wn 



n'=0 



(20) 



Using ([5]) and pTI) , chains are cloned and pruned if their 
weight is above 3Z„ and below Z„/3, respectively. 

Results of \\iZN{q,h) + jiooN plotted against N are 
shown in Fig.SJJa) for various values of b. For small h the 
curves are close to the curve for h — \. As 6 increases, 
the initial (small- A^) parts of these curves are straight 
lines with less and less negative slopes. In this regime the 
polymer is stretched. As long as these slopes are negative, 
the straight lines will intersect the curve for 6 = 1 at some 
finite value of A^, say Nc{b), i.e. for the finite system of 
chain length Nc{b) the corresponding effective transition 
point is b. For A^ > Nc{b), the values of In ZN{q, b)+iJ.ocN 
must deviate from the straight lines { see Refs. [2l| and 
[31| } for the detailed explanations. Since the curve for 
& = 1 becomes horizontal for A^ — >■ oo, the true phase 
transition occurs at that value of b for which the straight 
line in Fig. IH^a) is also horizontal. This can be estimated 
very easily and with high precision, giving for q = 1.5 
our final estimate be « 1.856(1). 

To clarify that the transition is indeed a first-order 
phase transition, one can study the histograms of x and 
m since PERM gives direct estimates of the partition sum 
and of the properly normalized histograms. The general 
formula of the histogram is 



P,,b{m,x)= J2 9" '" 

walks 



OfYijTn'^XjX' 



(21) 



Reweighting histograms obtained with runs performed 
nominally at q' and b' is trivially done by 



P,,l,{m,x) = Pg.^b'{m,x){q/q')"\b/b'y 



(22) 



Combining results from different runs can then be either 
done by selecting for each {m, x) just the run which pro- 
duced the least noisy data (which was done here in most 
cases), or by assuming that the statistical weights of dif- 
ferent runs are proportional to the number of "tours" [14!| 
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FIG. 4. (a) In .^jv(g, b) + fiooN for d = 3,q = 1.5, and for vari- 
ous values of b. The value fioo = — 1.7530±0.0003 used in this 
plot was obtained from dense limit simulations on finite lat- 
tices [3l|]. (b) Histograms of the end point distance P{x) ver- 
sus x/N for q = 1.5. Biases were adjusted so that both peaks 
have equal height: b = 1.4040 (iV = 500), 1.4925 (iV = 1000), 
1.5386 {N = 1500), 1.5658 {N = 2000), 1.5855 (iV = 2500). 
Normalization is arbitrary. The peak at x/N « corresponds 
to the collapsed phase, the other one to the stretched phase. 
Adapted from Ref. [s^- 



which contributed to Pq^i,{m,x). Note that for conven- 
tional Metropolis-type Monte Carlo algorithms, it is not 
trivial to combine MC results from different temperatures 



32|. 



since the absolute normalization is unknown 

An example of histograms P{x) for fixed q = 1.5 and 6, 
plotted against x/N are shown in Fig. Hljb) for A^ = 500, 
1000, 1500, 2000, and 2500. The value of b is determined 
such that the two peaks have the same height for each 
N, i.e., bc{N) = b is the effective transition point for 
the finite system of size A^. In addition the normaliza- 
tion factor is chosen arbitrarily to make all peaks having 
similar height for convenience. Using (P^. each curve in 
Fig. |4ljb) contributed by the properly reweighting data 
from different runs for various values of b. Obviously, 
with increasing A^, we see that the distance between 
two peaks increases and the minimum between the peaks 
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FIG. 5. Rescaled mean square end-to-end distance 

{R^)/{2£i,N^'^) plotted against the chain length A'^ for semi- 
flexible chains with It — 1 and various values of qi, on a log- 
log scale. Here u ~ 0.588 is the Flory exponent for SAWs in 
d = 3. Adapted from Ref. [H]. 



shrinks to zero. This gives a strong evidence for the first- 
order transition. Notice that a double peak structure 
with decreasing minimum alone would not be a conclu- 
sive proof, as shown e.g. by the 0-point in dimensions 
d > 4 [SjI and by some non-standard percolation mod- 
els fH. 



VI. SEMIFLEXIBLE POLYMER CHAINS 

Based on a Flory-like treatment [l|, [3a], for a chain 
with n units of the Kuhn length ix, and diameter d 
randomly linked together such that the contour length 
L = Nil, = n£x (there are {N + 1) monomers in the 
chain and connected by the bond length £},), the effec- 
tive free energy of such a semiflexible chain contains two 
terms as follows. 



A^ « Rl/ilKL) + v,Rl [{L/eK)/Riy 



(23) 



The first term is the elastic energy which is obtained 
by treating the chain as a free Gaussian chain, hence 
one can immediately write down the probability of the 
end-to-end distance Re which agrees with the Gaussian 
distribution. Therefore, the elastic energy is simply the 
logarithm of this distribution. The second term is the 
repulsive energy due to pairwise contacts where the sec- 
ond virial coefficient V2 = i^d, the density of monomers 
p = n/Rl = LRI/£k and the volume V = R^. Minimiz- 
ing AF with respect to Re, one obtains the Flory- type 
result for self- avoiding walks as L — ;> oo {N — > oo) 



Rp 



{v2/eK)'^'L'/'' = {iKdy/'{Ni,f/'' . (24) 



The minimum contour length L where the exclusive vol- 
ume is effective, i.e. the second term in (|23p is negligible 
in comparison with the first one if A^ < A^* , and using the 



scaling law of the square for the end-to-end distance of 
a Gaussian chain, i?g = ij^L = ii^iijN, the upper bound 
of the chain length for describing the Gaussian chain is 
obtained with 



N* = i%/{ibd^) 



(25) 



As L < £k, the chain shows a rod- like behavior, the lower 
bound of the chain length for the Gaussian chain is given 
by ik/ib- Therefore, the intermediate Gaussian behavior 
should only exist for 



4/4 <N <N* 



(26) 



For a linear semiflexible polymer chain (d = £b) under 
good solvent conditions, one would expect to observe 
both a crossover from rigid rod-like behavior to almost 
Gaussian random coils, then a crossover to self- avoiding 
walks when the chain stiffness varies. 

In order to verify the prediction, it requires an effi- 
cient algorithm to generate sufficient samples for very 
long semiflexible chains since the results should cover the 
linear length scales in the three different regimes. PERM 
was first applied to this in [3a]. The model described 
below had indeed been studied by means of PERM al- 
ready in [37,], where however most emphasis was put on 
the question whether the collapse transition changes from 
second to first order as the stiffness is increased. This was 
predicted by mean field theories [38] . The simulations 
in [33] supported the prediction, but were dangerously 
close to the significance limit discussed in Sect. 2. 

The above scaling relations for chains without self at- 
traction were studied in [361 ] . Semiflexible polymers were 
there modelled by SAWs on the simple cubic lattice, with 
a bending energy efc(l — cos0). Here is the angle be- 
tween the new and the previous bonds (only 9 = and 
6 = ±7r/2 are possible on a simple cubic lattice). The 
partition function of the SAWs of A^ steps with A^bcnd 
local bends (where 9 = ±7i'/2) is 



^A^.A^bond [lb 



ilb 



E GiN,N^ 

config. 



ad)gb 



A^b 



(27) 



where qb = exp(/3eti) is the appropriate Boltzmann fac- 
tor {qb = I for ordinary SAWs), and C(Af, A^bcnd) is the 
total number of chain configurations containing (A^ -1-1) 
monomers and A'bcnd local bends. 

In the simulation, the walk of length n — 1 at the n^^ 
step can be guided to either walk straight ahead in any 
direction, or make an _L-turn. Of course, it is only al- 
lowed to walk to the free nearest neighbor sites of the 
^th nionomer. The ratio of probabilities between the for- 
mer case and the latter case is chosen as l/g&. Since 
the stiffness of the chain is controlled by qb, we give 
less probability to make an L-turn as qb becomes smaller 
which corresponds to the case that the chain is stiffer. 
Results of the rescaled mean square end-to-end distance 
{RD / {2(.bN'^'') plotted against the chain length A^ up to 
A^ = 50000 for 0.005 < qb < 1.0 are shown in Fig. For 
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z^D+l 



z = 



FIG. 6. Schematic drawing of a polymer chain confined be- 
tween two waUs located at z = and z = D + 1. For our sim- 
ulations, chains are grown from the starting point {xo, yo,zo)- 
Here xq and yo are fixed but zq — 1,2, . . . ,D. 



stifFer chains, namely for smaller values of qt, we do see a 
rod-like regime at the beginning for not very long chains 
then a cross-over to a Gaussian regime, and then finally 
the excluded volume effect becomes more important for 
very long chains, and a horizontal plateau is developed. 
For very small qb, although the maximum chain length is 
up to 50000, it does still not yet reach the SAW regime. 
However, this is the first time that one can give evidence 
for the existence of the intermediate Gaussian coil regime 
(|26p by using computer simulations. 



VII. POLYMERS IN CONFINING 
GEOMETRIES 



Polymers Confined between Two Parallel Hard 
Walls 



It is a challenge to verify the theoretical scaling pre- 
dictions for single polymer chains of length N confined 
between two parallel hard walls with distance D away 
from each other (Fig. 16]) due to the difficulty of producing 
long polymer chains by MC simulations and the existence 
of very large finite-size corrections. For unconstrained 
SAWs, it is well know that the asymptotic scaling be- 
havior is reached rather slowly with correction terms de- 
creasing only as N~°-^ |39l - l4l| . Therefore, in addition to 
SAWs, we studied also the Domb- Joyce (DJ) model [42[ 
with V = 0.6 (where convergence to asymptotia is much 
faster [40,|4l|). 

In the DJ model, polymers are described by lattice 
walks where monomers sit at sites, connected by bonds 
of length one, and multiple visits to the same site are 
allowed, (i.e. the polymer chain is allowed to cross itself), 
but the weight is punished by a repulsive energy e > 
for any pair of monomers occupying the same site. Each 
pair contributes a Boltzmann factor v = exp(— /3e) to the 
partition sum. Thus, the partition sum of a linear chain 
consisting of A^ + 1 monomers is given by 



configs. 



(28) 



where the sum extends over all random walk (RW) con- 
figurations with N steps, < ti < 1, and m is the to- 
tal number of monomer pairs occupying a common site, 
m — X]i<7" '^xiXj (xj denotes the position of the monomer 
i). For u = 1, it corresponds to the ordinary RW. For 
f = it is just the SAW model. In the thermodynamic 
limit where A — > cx), the DJ model is in the same uni- 
versality class of SAW for all u < 1. There is a "magic" 
value of the interaction strength v — v* ^ 0.6 where cor- 
rections to scalin g ar e minimal and asymptotic scaling 
is reached fastest [40, |41| . In the renormalization group 
language, the flow speed of the effective Hamiltonian ap- 
proaching its fixed point depends on v. Moreover, it is 
approached from opposite sides when v < v* and when 
V > V* , with V* = 0.6. 

There exist important theoretical predictions for the 
monomer density profile p{z) and the end monomer den- 
sity profile Pe{z) near the wall given by Eisenriegler et 
al. Ei, [H as follows: 



P(^) 



.V-^a 



and 



P.ndiz) ^ Z^^--^' '^1- ^ Z 



^0.814(6) 



(29) 



(30) 



where z is the distance from the wall and 7*-^-' is the 
entropic exponent for 3D SAW with one end grafted on 
an impenetrable wall. One should also expect that the 
density near the walls is proportional to the force per 
monomer /. Indeed it was shown by Eisenriegler [46[ 
that 



lim/fc4fl 



= B 



f 



knT 



^B- 



V'Ap.o 



D 



-1-1/1/3 



(31) 



with B being a universal amplitude ratio. For ideal 
chains one has B — 2, while for chains with excluded 
volume in 4 — e dimensions one has _B w 2(1 — 6ie) with 
61 — 0.075 [43]. In three dimensions this gives the pre- 
diction B ss 1.85. 

In order to check the above mentioned theoretical pre- 
dictions, we simulate the SAW model and the DJ model 
on the simple cubic lattice with the confinement of a slab 
with width D by using PERM with 6-step Markovian an- 
ticipation. For estimating the monomer density profiles 
p{z) we only count those monomers in the central part 
of the chain, excluding 10% on either side to avoid er- 
rors from the fact that (P5| should hold only far away 
from the chain ends, for monomer indices n satisfying 
D^ <^n<^N~D^ (we should mention that N/D"^ > 10 
for all data sets). Results of p{z) obtained from the sim- 
ulations are normalized such that X]z=i P(-^) ~ ^- Since 
we simulate single polymer chains between two walls at 
z = and z ^ D + 1, we can assume that 

(32) 
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FIG. 7. Results of the monomer density profiles p[z) obtained 
for the DJ model, (a) Rescaled values of monomer density 
(D + V)p{z) plotted against ^ = ^/{D + 1). The function 
/o(C) = 18.74(^(1 - C))^''"^- (b) The same data as in (a), 
divided by foiX), plotted against a modified scaling variable, 
^i = (2 + 5)/(D + 1 + 2(5) with 5 = 0.04. The prefactor in ^ 
for a = and a = (D + 1) — >■ oo is 0.71(3). Adapted from 
Ref. O. 



where the constant A = 18.74 is determined by normal- 
ization. We plot the rescaled values of the monomer den- 
sity {D + ^)p{z) against £, in Fig. El^a) for the DJ model. 
It looks like that the scaling law ([2^ is satisfied and our 
data are described by the function /o(C) quite well for 
z & [0, D -I- 1]. But, we actually miss the important in- 
formation near the two walls in such a plot. A prefactor 
on the right hand side of (j31]) is probably not a constant. 
In order to give a precise estimate of the amplitude B 
([5T|) we introduce here an "extrapolation length" 5 as 
suggested in [4^ so that the scaling variable ^ is replaced 

by 



6 



z + 5 



D + l + 25 



(33) 



Using the same data of p{z) but divided by 
fo{£,s), the best data collapse shown in Fig. Eljb) 
is obtained by taking S = 0.04. It leads to 



\im^^o,D^cc> D^+^/'"> z-^/'"> p{z)/A = 0.71(3). Since the 
extrapolation length 6 = 0.04 for the DJ model is much 
smaller than ^ = 0.15 for SAWs {see Fig. 7 in Ref. (H), 
it gives a first indication that corrections to scaling are 
indeed smaller in the DJ model. Using ((3T|) . it gives 
B = 1.70 ± 0.08. This is only 2 standard deviations 
away from the renormalization group expansion predic- 
tion or Ec = 4 — d expansion prediction B — 1.85 of 
Eisenriegler [46|, which we consider as good agreement. 



B. Escape Transition of a Polymer Chain from a 
Nanotube 



The confinement or escape problem of polymer chains 
in cylindrical tubes of finite length has the merit that it is 
potentially very relevant to experiments and applications 
such as the problem of polymer translocation through 
pores in membranes and the study of DNA confined in 
artificial nanochannels [a-[g|. The following treatment is 
based on [H-ilj. 

Considering a polymer chain of length N with one end 
grafted to the inner wall of a cylindrical nanotube with 
finite length L and diameter D under good solvent condi- 
tions, the chain configuration is compressed uniformly as 
D decreases or N increases, but L is fixed. Beyond a cer- 
tain compression force, the chain configuration changes 
abruptly from a homogeneously stretched and confined 
state (imprisoned state) to an inhomogeneous state (es- 
caped state) where polymer chains form a flower-like con- 
figuration with one stem confined in the tube and a coiled 
crown outside the tube (see Fig,[S]). This abrupt change 
implies a first order transition. Since the theory based 
on the blob picture failed to predict the transition from a 
homogeneous state to an inhomogeneous state, the Lan- 
dau theory approach is used for describing such a first 
order transition including the metastable states. In the 
Landau theory approach, all configurations are subdi- 
vided into subsets associated with a given value of an 
appropriately chosen order parameter s that allows to 
distinguish between different states or phases. The full 
partition function of the system is therefore obtained by 
integrating over the order parameter: 



Z = 



exp(-F) = / dsexp[-$(s)] , 



(34) 



where $(s) is the free energy of a given set, and is there- 
fore a function of the order parameter. Here the order 
parameter s is defined by the stretching degree, i.e. the 
ratio between the end-to-end distance of monomer seg- 
ments which are still confined in the tube, i?imp, and 
the number of monomers confined in the tube, iVjmp- As 
shown in Fig. ^ we see that near the transition point the 
Landau free energy function has two minima, the lower 
minimum is associated with the thermodynamically sta- 
ble state, which corresponds to the equilibrium free en- 
ergy (either Fimp or i^csc) of the system, while the other 
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FIG. 8. Schematic drawings of a flexible polymer chain of length A'^ grafted to the inner wall of a tube of length L and diameter 
D at the transition point, (a) As the chain is fully confined in the tube (in an imprisoned state), it forms a sequence of 
Ub — ND~^''^ blobs in a cigar-like shape, here // = j/a is the Flory exponent in d = 3. (b) As one part of the chain escapes from 
the tube (in an escaped state), it forms a flower- like configuration which consists of a "stem" containing A'tr monomers and a 
"crown" containing A'^ — A'tr monomers. 
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FIG. 9. The Landau free energy per monomer, 

$(7V, L, D, s)/N plotted against the order parameter s near 
the transition point for the tube of length L = 1600 and di- 
ameter D = 17. 



minimum corresponds to the metastable state jSOl ISlf . 
At the transition point, both minima are of equal depth. 

In our simulations, we describe the grafted single poly- 
mer chain confined in a tube by SAWs of N steps on a 
simple cubic lattice with cylindrical confinement {0 < 



X < L,y^ 



D^/A}, and the first monomer is at- 



tached to the center of the inner wall of the tube. Taking 
the advantage of PERM that the associated weight of 
each generated configuration is exactly known, we intro- 
duce a new strategy in order to obtain sufiicient sam- 
plings of the flower-like configurations in the phase space 
as follows: We first apply a constant force along the tube 
to pull the free end of a grafted chain outward to the 
open end of the tube as long as the chain is still con- 
fined in a tube, and release the chain once one part of 
monomer segments of it is outside the tube. Varying 
the strength of the force, we obtain flower-like configura- 
tions containing stems with various stretching degree of 
monomer segments which are still confined in a tube if 
the length N is long enough. The contributions for the es- 
caped states are therefore given by properly reweighting 
these configurations to the situation where no extra force 
is applied. This is done by using biased SAWs (BSAWs) 



on a simple cubic lattice with finite cylindrical geome- 
try confinement, similar to the model in (|16p . but we use 
here q = 1 to describe the good solvent condition. 

With PERM, the total weight of a BSAW of N steps 
(A^ + 1 monomers) is Wb{N, L, D) — H^^qWu with w„ = 
&(^"+i-=""VPn for n > 1 and wo{N,L,D) = 1. p„ is 
chosen as in (|18l) . The estimate of the partition sum is 
given by 



Zb{N,L,D) 



- E 



Mb 



Wb{Cb) 



(35) 



configs.t^Cb 



where a set of configurations is denoted by Cb- Thus, 
each configuration of BSAWs with the stretching factor 
bu contributes a weight l^C^) [N, L, D) for a BSAW of N 
steps with b = 1 confined in a finite tube of length L and 
diameter D: 

W^^Hn LD)^l Wb,{N,L,D)/bl'^^'-^' , XM<L 
yi^,i.,u) \^^^^^^L,D)/bi , XN>L ' 

(36) 
where index k labels runs with different values of the 
stretching factor b. Combining data runs with different 
values of 6, the final estimate of the partition sum is 

Z{N^L,D)^^Y. E W^'^^N^L^D) (37) 



k configs.^Cb, 



here M is the total number of trial configurations. 

The distribution of the order parameter, 
P{N,L,D,s) ex H{N,L,D,s), is obtained by accu- 
mulating the histograms H{N,L,D,s) of s, where 
H{N, L, D, s) is given by. 



H{N,L,D,s) - iTJ:kH^'HN,L,D 



\I l^k Z^ 



M 



configs 



eCb,W^''HN,L,D,s')Ss 



and the partition sum of polymer chains confined in a 
finite tube can be written as 



ZiN,L,D) = J2H{N,L,D,s) 



(38) 



in accordance with (P7)) . Thus, one can also double check 
the results of the partition sum. 
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The Landau free energy ^{N, L, D, s) here is the ex- 
cess free energy related to the polymer chains with 
one end tethered to an impenetrable flat surface, i.e. 
cP{N,L,D,s)^ ~\n[P{N,L,D,s)/Zi{N)] {Z,{N) ~ 
^Nj^ji-i f^) Results shown in Fig.[9]are for L = 1600, 
D = 17, and for N/L = 5.5, 5.7, 5.9. This shows that 
the information about metastable states can also be ex- 
tracted from the simulations with PERM. 



VIII. PERM FOR BRANCHED POLYMERS 
WITH FIXED TREE TOPOLOGIES 



g{n), for n, > 0, 



lin) 



(n + 4.0)/(n + 1.3) 
(n + 0.6)/(n + 3.9) 



outward direction 
otherwise 



However, the strength of this bias can be adjusted 
by trial and error. 

• The population control (pruning/cloning) is done 
in the same way as explained in Sect. |TT] that at 
the step n, two thresholds W^ and W~ are pro- 
portional to the current estimate weight Z„, e.g., 
W+ = 2,Zn and W' = 0.5Z„. 



In this section we shall discuss two types of branched 
tree-like polymers: Star polymers (where all branches 
emanate from one single point) and "bottle brushes" 
where side chains of common lengths are attached to a 
backbone at regularly spaced points. 

To be concrete, let us consider the simplest case of 
a branched polymer, a star polymer where / arms are 
grafted to a single branch point, and all arms have the 
same length N . 

As a linear chain is built by using PERM, at each step 
one monomer is added to the built chain until the chain 
has reached its maximum length N or it has been killed 
in between. For growing a star polymer we have to be 
aware that not only the interactions between monomers 
in the same arm have to be considered but also the in- 
teractions between monomers on different arms have to 
be taken into account. If one arm is grown entirely be- 
fore the next arm is started, it will lead to a completely 
"wrong" direction of generating the configurations of a 
star. However, it is straightforward to modify the basic 
PERM algorithm such that all / arms of a star poly- 
mer are grown simultaneously [43l . |53|. The multi-arm 
method is explained as follows: 

• A star polymer is grown from its branching point 
(center). 

• / growth sites {xi,...x^} are considered at the 
same time. A monomer is added to each arm step 
by step until all arms have the same length, then 
the next round of monomers is added. As all the 
monomers in a star are numbered, it is similar as 
growing one linear chain from the 1*** monomer to 
the N^^^ monomer (see Fig. ^. iVmax ^ Nf + 1 
if the center is singly occupied or A^max = N f + / 
if the center is /-folded occupied. 

• A bias is given to guide the growth of arms into 
outward direction with higher probability. The 
strength of this bias is adjusted in the way that 
it increases with / but decreases as the length of 
arms becomes longer since there is more space in 
a dilute solvent for adding the next monomer. For 
example, we can choose the bias as a function of n. 



A. Star Polymers 

For single star polymers composed of / arms of length 
N each in a good solvent, the partition sum and the rms 
center-to-end distance scale as follows: 



and 



7(1) 



Rnj 



-fNjsf-rf- 



AfN 



2u 



(40) 



(41) 



where the critical fugacity fi,^ and the Flory exponent 
v are the same for all topologies but the entropic expo- 
nent jf depends on each topology [54| . In two dimen- 
sions, 7/ can be calculated exactly by using conformal 
invariance ^54], but there are no exact results for the /- 
dependent power law for jf, and also not for the swelling 
factor Af. Therefore, computer simulations are needed 
for a deep understanding of star polymers. Due to the 
difficulty of simulating the star polymers with many arms 
f and of long arm length N by both MC simulations (Sa - 
l60J and molecular dynamics [6ll, |62| , and because of the 
lack of precise estimates of the exponents given in (j40|) 
and (|41l) . PERM with multi-arm growth method as ex- 
plained above was developed |43|. With this algorithm, 
high statistics simulations are obtained for star polymer 
with arm number up to / = 80 and arm length up to 
N = 4000 for small values of /. 

For our simulations of single star polymers in a good 
solvent, we use the Domb- Joyce model with the interac- 
tion strength v* = 0.6 on the simple cubic lattice (see 
Sect. IVIIAp . It allows us to attach a larger number of 
arms to a point-like center of stars, and thus additional 
considerations of the corrections to scaling terms when 
a finite size core is used are avoided. Two variants for 
studying star polymers are used in our simulations. In 
one variant the center is occupied by one monomer, and 
in the other variant the center is occupied by / monomers 
as shown in Fig. 1101 Since the partition sum is estimated 
directly by PERM, the exponents 7/ can also be deter- 
mined easily according to (j40|). 

In Fig. 111! we present results of 7f from our simulations 
and from previous studies |56l . 1571 |63| for comparison. 
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(a) 
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arm 2 



(b) 




arm 1 
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FIG. 10. Schematic drawings of a star polymer consisting of three arms (/ = 3) of length A*' = 3 each. The center is singly 
occupied in (a) and /-folded occupied in (b). Those numbers show the order of monomers which is added into the star polymer 
by using a chain growth algorithm. 
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FIG. 11. Exponents 7/ plotted against /. The full line is 
just a polygon connecting the points, and the dashed line 
is a fit with the large-/ behavior as predicted by the cone 
approximation (|42[) . Results obtained in Ref. [Sa. ISTI. |63|| are 
shown for comparison. In the inset, we show those results for 
small /. Adapted from Ref. [43 |. 



The theoretical prediction for the scaling 
large / by the cone approximation [59l l64 



7/ - 1 -- / 



-3/2 



law of 7/ for 
is 

(42) 



For small /, our results are in good agreement with the 
previous studies. For large / the best fit with a power law 
7/ — 1 ~ — (/ — 1.5)^ would be obtained with z w 1.68, 
which is not too far off the theoretical prediction (j42|) but 
the prediction is also not exact. 

After we have obtained quite reliable estimates of ^ao , 
v, and 7/ for single star polymer in a good solvent, we 
extend our study to a more complicated s yste m where 
two star polymers interact with each other [53| by using 
the same model and the same algorithm. It is well under- 
stood that interactions between both linear and branched 
polymers are soft in the sense that they can penetrate 
each other and the effective potential is a rather smooth 
function of their distance. For star polymers, there are 
some contradictions between results in the literatures. Is 
the potential between two central monomers at large dis- 
tance a Gaussian potential or has it a Yukawa tail? Since 



we were able to simulate star polymers up to / = 80 
arms, we expected that we would give a clear answer. 
This was the main motivation to study the effective po- 
tential between two star polymers [SSj . 

Witten and Pincus [6J| point out that the scaling 
of the partition sum of a star with / arms and arm 
length N each (|l(I)l . together with the assumption that 
(1) 



4'/(^)/ 

fixed /, i.e 



Z 



N,f 



is a function of a; = f/^g only for any 



4'?/(0 






^fir/Rg) 



implies that 



V(r) « Fwp(r) = bfln{afRg/r) 



(43) 



(44) 



where r is the distance between the two central 
monomers, and 

bf = {2jf - 72/ - l)/!^ for 1< r < i?g . (45) 

According to our results shown in Fig. Ill[ instead of the 
scaling &/ ~ f^^^, a power law gives 6/ w 0.27/^'^®. 
However, both af and bf should be universal and should 
not depend on the specific microscopic realization. 

There are two methods for estimating Z^'^'{r) in our 
simulations; 

(a) Two independent star polymers are grown simulta- 
neously, and Z^'^'(r) is computed by counting their 
overlaps at different distance r. Here Z'-^^r) and 
Z^^' (r) are estimated in the same run, which gives 
rather accurate results for the potential V{r) for 
very large distances r and large N. For small dis- 
tances r, the ratio Z^'^\r)/ [Z^^)(r)] would be in- 
distinguishable from zero. 

(b) Two star polymers are grown at fixed distance r 
with the mutual interactions taken into account 
during the growth. This allows us to measure 
2'(2)(j^) down to very small distances r and large 
A'^. For large distances r, it gives very bad results 
of the potential V{r) since it is obtained by sub- 
tracting the (nearly equal) free energies obtained 
in two different runs. 
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FIG. 12. (a) The effective potential V{r) for / = 18, plotted 
against Rg in a semi-log scale. The solid curve shows (|44|l . 
and the dotted curve is a Gaussian, (b) Rescaled radial Mayer 
functions against r/Rg for several values of /. Curves are 
obtained from (|48|l . with fitted parameters a/, c/, df and r/. 
Adapted from Ref. [s^ l 



For the data analysis, we use the data either from the first 
method or the second method, or use the combination 
from both. 

We present the effective potential V{r) between two 
star polymers of / = 18 arms in Fig. [T^a). For r <C 
Rg, V{r) follows the prediction given in (|33]), which is 

For r :^ Rg the MC data 



shown by the solid curve 

can be approximated by a parabola, i.e. V{r) is roughly 

Gaussian 



V{r) « FGauss(r) = c/e'*^'-'/^' 



(46) 



Here we conjecture that Cf and df are universal. In or- 
der to describe the effective potential V{r) for the whole 
region of r, we propose that 



V{r) = — In 



where ry is an additional parameter for every /, and 

V{r) > for all r. As r ^ oo, V{r) = VGauss(V)[l + 



0{r-^f)] 66]), while Vir) = Vwp(r-)[1 + ©(r^)] gl) as 
r — > 0. In figHHb), we plotted the rescaled radial Mayer 
function, 



ir/Rl)'fMir) = {r/Rg)\a-eM-Vir)]) 



(48) 



against the rescaled distance r/Rg. Our results are in 
good agreement with the simulations of [g^ but do not 
agree with the results in J66|. 



B. Bottle-brush Polymers 

The so-called bottle-brush polymer consists of one long 
molecule serving as a backbone on which many side 
chains are densely grafted. As the grafting density a in- 
creases, the persistence length of the backbone increases. 
The bottle-brush polymer has the form of a rather stiff 
cylindrical-like object. If the backbone is very short but 
side chains are very long, it should behave like a star poly- 
mer. If the backbone is very long, the structure becomes 
more complicated. One would expect that those side 
chains in the interior of the bottle-brush are all stretched 
and show the same behavior, but those at the two ends 
behave as a star. In order to understand the structure 
of bottle-brush polymers and check the scaling behavior 
of long side chains in comparison with theoretical predic- 
tions [67|, we focus here the bottle-brush polymers of a 
rigid backbone and flexible side chains. 

For our simulations, we use a simple coarse-grained 
model. The backbone is treated as a completely rigid 
rod, and side chains are described by SAWs with near- 
est neighbor non-bonded attractive interactions between 
the same type of monomers and repulsive interactions 
between the different type of monomers. A general for- 
mula for the partition function for bottle-brush polymers 
consisting of one or two chemically different monomers is 
therefore given by 



Z 



1 Iab 



(49) 



where q = exp(— /?e) (we assume that the attractive inter- 
action eAA = ^BB = e), qAB = exp(-;3eAB) {^ab is the 
repulsive interaction between monomer A and monomer 
-B), and tuaa, ^ubbj it^ab are the numbers of non-bonded 
occupied nearest neighbor monomer pairs AA, BB and 
AB, respectively. For q — 1, all side chains behave 
as SAWs. For q < q^ it corresponds to the good sol- 
vent condition, where q^ = exp(— e/fcsTe) ~ 1.3087 at 



the point [IJ]. For q > q^, it corresponds to the 
poor solvent condition. As qAB = 0, it corresponds 
to a very strong repulsion between A and B, while for 
QAB — Q the chemical incompatibility vanishes {recall 
that XAB oc eab - {<^AA + esB)/2}. The grafting den- 
sity a is defined by ct = Uc/Nh where Uc is the number 
of side chains and Nb is the number of monomers in a 
backbone. Here only the results of bottle-brush poly- 
mers consisting of one kind of monomers under a very 



15 



(b) 





FIG. 13. (a) A schematic drawing of growing a bottle-brush polymer step by step, (b) A snapshot of the configurations of 
bottle-brush polymers consisting of A^i, = 128 backbone monomers, N = 2000 monomers in each side chain, and the grafting 
density a = 1/4 under a very good solvent condition generated by PERM. 



good solvent condition are presented in order to show the 
performance of the algorithm. Other applications can be 
found in [IMOl- 

We extend the algorithm for simulating star polymers 
to bottle-brush polymers. As shown in Fig. [TST a). a 
bottle-brush polymer is built by adding one monomer 
to each side chain at each step until all side chains have 
the same number of monomers. Then we start to add the 
second run of monomers, i.e, all side chains are grown si- 
multaneously. The bias of growing side chains was used 
by giving higher probabilities in the direction where there 
are more free next neighbor sites and in the outward di- 
rections perpendicular to the backbone, where the second 
part of bias decreases with the length of side chains and 
increases with the grafting density. A typical configu- 
ration of bottle-brush polymers consisting of Ni, — 128 
backbone monomers, N — 2000 side chain monomers, 
and with grafting density cr = 1/4 under a good solvent 
condition is shown in Fig. llSf b) where the total num- 
ber of monomers is iVtot = 128 + 2000 x 32 = 64128 
monomers. 

For checking the scaling law of side chains, we in- 
troduce the periodic boundary condition along the di- 
rection of the backbone (-l-z-direction) to avoid end ef- 
fects associated with a finite backbone length. The 
square of the average height of a bottle-brush polymer, 
Rl{N,a) = {Rl^{N,a) + Rly{N,a)) is estimated by 
taking the average of the mean square backbone-to-end 
distance in the radial direction for all side chains. In 
Fig. [IHa) we plot R\{N,a) divided by iV^" versus N 
for Nh — 32, 64, and 128, for various values of grafting 
densities a. The value of v is given by the best esti- 
mate for 3d SAW by PERM [4^. We see that those 
curves of the same grafting density a coincide with each 
other. Increasing the grafting density a, it enhances the 
stretching of side chains. As u — ^ 0, we should ex- 
pect a mushroom regime where no interaction between 
side chains appears. As a is very high, the scaling pre- 



diction obtained by extending the Daoud-Cotton [7l| 
"blob picture" |72h75| from star polymers to bottle-brush 
polymers is Rl{N,a) ex (j'^{^-^)/{^+'^)n^'^/{i+^') , Thus, 
we can give the cross-over scaling ansatz as follows for 

N ^ CO, 



Ri{N,a)^N"^R'{^) 



with 



R^iri) 



2(l-i.)/(l+7.) 





oo 



(50) 



(51) 



where 77 = aN^ . 

After removing those unphysical data due to the arti- 
fact of using periodic boundary condition in the regime 
where Rh{N,a) > Nb/2, we plot the same data of 
R\{N , a) / N'^'^ but rescaled the x-axis from N to rj — 
aN" according to the scaling law (15(11) . We see the nice 
data collapse in Fig. [T4lb). In this log- log plot, the 
straight line gives the asymptotic behaviors of the scaling 
prediction (j50p for very large rj. As 77 increases, we see 
a cross-over from a 3D SAWs to a stretched side chain 
regime but only rather weak stretching of side chains is 
realized, which is different from the scaling prediction. 
However, this is the first time one can see the cross- 
over behavior by computer simulations. This cross-over 
regime is far from reachable by experiments. 



IX. PERM WITH CLUSTER GROWTH 
METHOD 

It is generally believed that lattice animals, lattice 
trees, and subcritical percolation are good models for 
studying randomly branched polymers and they are in 
the same universality class. There exist several efficient 
algorithms, e.g., Leath algorithm [7g, Swendsen-Wang 
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FIG. 14. (a) Log-log plot of rescaled mean square height 
Rl{N, a)IN'^" versus iV (a) and r] = aN" (b) with v ^ 0.588. 
Results are obtained for three choices of Nt and several 
choices of the grafting density a as indicated. Those unphys- 
ical data {Rh > 0.5Nb) due to the artifact of using periodic 
boundary condition are removed. The slope of the straight 
line corresponds to the scaling prediction. Adapted from 
Ref. M. 



algorithm ^], etc. for studying the growth of percola- 
tion clusters near the critical point, but they all become 
inefficient far below it, because the chance for growing a 
large subcritical cluster by a straightforward algorithm 
decreases rapidly with N. Obviously we need some sort 
of cloning, and since this will probably lead also to fluc- 
tuating weights, one might need some pruning. 

Cloning and pruning needs first some estimate for the 
weight of a cluster that is still growing. Moreover, it will 
turn out that growing clusters can have, depending on 
their detailed configurations, very different probabilities 
to grow further. Thus, in addition to a weight we might 
to need also a "fitness" that should depend on the weight 
but is not entirely determined by it. 

In the following discussion the algorithm is explained 
by considering the relationship between the site percola- 
tion and site lattice animals |78j . 

In any cluster growth algorithm |76| , a finished cluster 



with N sites and b boundary sites on a lattice is generated 
with probability 



P, 



Nb 



P''{l-P)\ 



(52) 



if each lattice site is occupied with the probability p. By 
definition of lattice animals all the clusters of same size 
N carry the same weight. Since the obtained percolation 
cluster is biased by the probability P/vb, its contribution 
to the animal ensemble is corrected by a factor 1/P/vh. 
Taking an average over the percolation ensemble, the par- 
tition sum of lattice animals consisting of N sites is given 
by 



Z 



N 



1 



{-—)=p-^^{l-p)-^) 



(53) 



As shown in Fig. I15[ now we consider a cluster with 
N sites, g growth sites and b boundary sites. At each 
of the growth sites the cluster can either grow further, 
or it can stop growing with the probability 1 — p. Thus, 
this still growing cluster gives a weight to a percolation 
cluster with N sites and (b + g) boundary sites as p^(l — 
p)''+s/ [p^(l — p)''] = (1 — p)^ ■ Taking an average over 
all clusters, we have 



Z 



N 



(l-p)fl 



>^(l-p)''+9 



)^p-''{il-p)-') 



(54) 



This is the same formula as given by (j53|) . but note that 
now we have included also those clusters which are still 
growing. 

Let us first point out this new variant of PERM: 

• The percolation cluster growth algorithm with stor- 
ing the growth sites into a queue in a first-in first- 
out list (the scheme of breadth-first) is used. 

• The population control is done by introducing a 
fitness function 



f^^Wn/il-pr^p-"il-p) 



-b-ag 



(55) 



with a parameter a to be determined empirically, 
and used 



/„>C+(/„), fn<C-{f„) 

as criteria for cloning and pruning. 



(56) 



• The depth-first implementation in PERM is still 
used here. Namely, at each time one deals with only 
a single configuration of a cluster until a cluster has 
been grown either to the end of the maximum size 
N or has been killed in between, and handles the 
copies by recursion. 

• The optimal value of the probability p is p < p^, 
and p ^ pc as N ^- oo. 
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FIG. 15. A still growing cluster with N — 7 sites, 6 = 6 boundary sites and g = 6 growth sites on a square lattice. 
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FIG. 16. Growing clusters generated in the (a) depth-first and (b) breadth-first implementations. In both cases, p — Pc ~ 0.5927 
and A*' = 4000. Occupied sites and growth sites are depicted by small red points and big black points, respectively. Adapted 
from Ref. [zi. 



This algorithm was developed more or less by trial and 
error, guided by the following considerations: 

We first test the two common ways for growing the 
percolation clusters, (a) Depth-first: growth sites are 
written into a first- in last out list (a stack), (b) Breadth- 
first: growth sites are written into a first-in first-out list 
(a queue). In order to avoid the mix up with the depth 
implementation in PERM, we use stack and queue to dis- 
tinguish these two methods. Two typical 2-d clusters of 
size N = 4000 and at the critical point of percolation 
p = Pc = 0.5925, growing according to these two meth- 
ods are shown in Fig. [1^1 At first glance, one would 
expect that the cluster growing by storing growth sites 
in a stack might be more efficient than that the growth 
sites stored in a queue, because the number of growing 
sites was about 3 times larger than that for the latter 
case. But the truth is, after a few generations the de- 
scendents generated from the former case will die. On 
the other hand, the fluctuations in the number of growth 
sites are much bigger in the former case, the weights in 
([54| will also fluctuate much more, and we expect much 
worse behavior. This is indeed what we found numeri- 
cally: Results obtained when using a stack for the growth 
sites were dramatically worse than results obtained with 



a queue. 

Second, we check whether the efficiency is affected by 
the chosen order of writing the neighbors of a growth 
site into the list. Studying the percolation cluster in two 
dimensions, one can use the preferences east-south-west- 
north, or east-west-north-south, or a different random 
sequence at every point. We found no big differences in 
efficiency. 

Third, it would be far from optimal to do the pop- 
ulation control as explained in Sect. |lll i.e. by us- 
ing two thresholds W^ on the current weights Wn = 
p^"(l — p)~^. This would strongly favor clusters with 
few growth sites, since they tend to have larger values 
of 6, for the same n, and have thus large weights. But 
such clusters would die soon, and would thus contribute 
little to the growth of much larger clusters. Therefore a 
proper fitness function /„ is needed. 

Finally, we have to decide the optimal values of p em- 
pirically. It is clear that we should not use p > Pc, be- 
cause it is subcritical percolation that is in the same uni- 
versality class of lattice animal. One might expect p <^ Pc 
to be optimal because only minimal reweighting is needed 
for small p. This is indeed true for small TV, but not for 
large N . In order to reach large TV, it is more important 
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that clusters grown with p <SiPc have to be cloned exces- 
sively, otherwise, they would die rapidly in view of their 
few growth sites. In Fig. [17] we present the errors of free 
energies Fjq = — In Z^v for various values oi p va d — 2 
and d = 8. The statistical errors always eventually de- 

1 /2 

crease as 1/ [CPU time] ' , hence we show there one stan- 
dard deviation multiplied by [CPU time] (measured in 
seconds), for different values of p. Thus, we can compare 
the accuracy between those runs on different computers. 
For d = 2 (Fig. [TWa')'). each simulation was done for 
-^max = 4000 (although we plotted some curves only up 
to smaller N , omitting data which might not have been 
converged). We see clearly that small values of p are good 
only for small N . As N increases, the best results were 
obtained for p ^ pc- The same behavior was observed 
also in all other dimensions, and also for animals on the 
bcc and fee lattices in 3 dimensions (data not shown). In 
Fig- fTWb'l. we see the analogous results for d = 8 and for 
-^max = 8000, showing the errors are much smaller than 
those in Fig. [TTT a') . Indeed, the errors decreased mono- 
tonically with d, being largest for d = 2. Using p slightly 
smaller than pc we can obtain easily very high statistics 
samples of animals with several thousand sites for dimen- 
sions > 2. Another quantity which can help to check the 
reliability of our data is the tour weight distribution (see 
Sect.|n|). In Fig. [T51 we show the two tour weight distri- 
butions for two-dimensional animals with 4000 sites, for 
p = 0.57 and for p — 0.47. We see that the simulation 
with p — 0.57 is distinctly on the safe side, while that 
for p = 0.47 is marginal. In the log-log plot, it is seen 
that the tail of the distribution P(ln>V) for p = 0.57 de- 
cays faster than l/W, thus the product >VP(lnyV) has 
its maximum where the distribution is well sampled. 

Error bars quoted in the following on raw data (par- 
tition sums, gyration radii, and average numbers of 
perimeter sites or bonds) are straightforwardly obtained 
single standard deviations. Their estimate is easy since 
clusters generated in different tours are independent, and 
therefore errors can be obtained from the fluctuations 
of the contributions of entire tours (notice that clusters 
within one tour are not independent, and estimating er- 
rors from their individual values would be wrong). 

In addition to site animals, this algorithm can also be 
applied to bond animals and lattice trees for studying 
randomly branched polymers. A bond animal is a clus- 
ter where bonds can be established between neighbor- 
ing sites (just as in SAWs), and connectivity is defined 
via these bonds: if there is no path between any two 
sites consisting entirely of established bonds, these sites 
are considered as not connected, even if they are near- 
est neighbors. Different configurations of bonds are con- 
sidered as different clusters, and clusters with the same 
number of bonds (irrespective of their number of sites) 
have the same weight [73| ■ Weakly embeddable trees are 
bond animals with tree topology, i.e. the set of weakly 
embeddable trees is a subset of bond animals, each with 
the same statistical weight. Strongly embeddable trees 
are, in contrast, the subset of site animals with tree-like 
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FIG. 17. Statistical errors of In Zn for lattice animals in d = 2 
(a) and d = 8 (b) for various values of p. To make the different 
runs comparable, errors are multiplied by the square root of 
the CPU time measured in seconds. The cluster size N is 
up to 4000 in (a) and up to 8000 in (b). The percolation 
thresholds are pc = 0.5927 in d = 2, and pc = 0.0752 in 
d = 8. Adapted from Ref. ^7^1. 
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FIG. 18. Log- log plot of distributions of tour weights P(lnVK) 
of 2d animals with iV = 4000, for p = 0.57 and p = 0.47, 
together with a straight line representing the function y = 
const I W . Adapted from Ref. [78| . 
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FIG. 19. (a) A site animal with 8 sites, (b) A site tree ("strongly embeddable tree"), (c) A bond animal which is not a tree, 
(d) A bond tree ("weakly embeddable tree"). 
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FIG. 20. Typical site lattice animals with N = 12000 on the square lattice in d = 2 (a), with A*' = 16000 on the bcc lattice in 
d = 3 (b). 



structure. All these definitions are illustrated in Fig. [T^ 



A. Non-interacting Lattice Animals in the Bulk 

The basic problem of lattice animals (site animals) is 
how to count the number of different animals of N sites 
precisely, i.e. the estimate of the corresponding partition 
sum. Two animals are considered as identical if they dif- 
fer just by a translation, but they considered as different 
if a rotation or reflection is needed to make them coin- 
cide. Two typical site animals consisting of A^ = 12000 
sites on the square lattice in d = 2 and with N = 16000 
sites on the body centered cubic (bcc) lattice in d = 3 
are shown in Fig. [201 

In the thermodynamic limit as A'^ — > oo, the number of 
animals (i.e . the microcanonical partition sum) should 
scale as [801 



Z 



N 



f.'^N'^'il + b.N-'^ + ■■■), 



and the gyration radius as 

Rn-^ N^a + bRN-^ + ---) 



(57) 



(58) 



Here fi is the growth constant (or inverse critical fugac- 
ity), and is not universal, while the Flory exponent v, the 
entropic exponent 9, and the correction exponent A [81| 
should be universal, bz and bn are non- universal ampli- 
tudes, and the dots stand for higher order terms in 1/A^. 
Results of the partition sum Z]\f and the mean square 
end-to-end distance _R^ for site lattice animals in d = 2 



are shown in Fig. [5TJ By taking the predicted value of 
6 = 1 and plotting IiiZn — aN + InA^ against A^, we 
should expect a curve which becomes horizontal for large 
N by adjusting values of a = In /i suitably. This is indeed 
seen for the central curve with error bar in the inset of 
Fig. [2lT a). but a precise estimate of /i is difficult because 
of corrections to scaling. Considering the first correction 
term in (|57l) and (|58|) . the correction exponent A, and the 
estimate of the growth constant /i and their error bars are 
all determined by the best straight line as A^~^ — ?> in 
Fig. [m Our estimate of a = In/x = 1.4018155(30) with 
A = 0.9(1) is ill perfect agreement with the exact enu- 
meration result [83] . The Flory exponent v is determined 
by the same way and our estimate v = 0.6412(5) is also 
in good agreement with the previous estimate by Monte 
Carlo simulations [83j . 

It is trivial to generalize the algorithm PERM with 
cluster growth method to lattice animals in higher di- 
mensions 2 < d < 9. Using the similar method of data 
analyses as shown in Fig. [21] for those results obtained 
in d = 3 to d = 7 (dc = 8 is the upper critical dimen- 
sion of lattice animals, where large corrections have to 
be taken into account). The relationship between the 
entropic exponent 9 and the Flory exponent v for the 
animal problem in d dimensions is predicted by using 
supersymmetry [84], 



{d-2)v + l 



(59) 



By plotting our data of the exponent v and {9 — l)/{d — 
2) against d in Fig. [521 we see that these two curves 
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FIG. 21. (a) Results of In Zjv - aiV + 6* In iV plotted against 
Af"^, and against A'^ (in the inset), and (b) results of R%/N'^" 
plotted against A^~ , and against A'^ (in the inset). The best 
estimates of a = In/i = 1.4018155(30), v = 0.6412(5) and 
A = 0.9(1) are given by the best straight lines. All data are 
for site lattice animals in d = 2. Adapted from Ref [781 ]. 
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FIG. 22. The critical exponents u and (6* — l)/(d — 2) against 
d. Adapted from Ref. [78|. 



for it, there exist no literature values to compare to our 
measurements. This is different for a second new expo- 
nent specific for the transition point, the cross-over expo- 
nent (j). If q is the Boltzmann factor for the monomer- wall 
interaction and qc is its critical value, then the scaling 
ansatz for the partition sum of a grafted animal near the 
adsorption transition is 



Z\^H,q)-^t^{qfN-'^n{q-q,)N^] 



(60) 



The most interesting prediction for was that is supe- 
runiversal, i.e. its value is independent of the dimension 
and (/) = 1/2 for all dimensions [85|- While this was veri- 
fied by the simulations for d = 3, 4 and 5, it was slightly 
violated (by 5 standard deviations) in d = 2 [7^. Ob- 
viously further investigations would be needed to settle 
this problem. 



coincide with each other. It shows that the Parisi-Sourlas 
prediction ([59]) is verified. 



Conformal Invariance and Animals Grafted to 
Wedges 



B. Lattice Animals Grafted to Surfaces 

For a-thermal walls (which represent only a geomet- 
ric barrier, without any other interactions) the leading 
behavior for TV — > cx) does not involve any new critical 
exponent [73|. This is no longer true, however, if the wall 
is attractive. In that case we expect a phase transition 
at a critical attractive energy beyond which the animal 
gets adsorbed to the surface, similar to the adsorption 
transition observed also for linear polymers ^52i|. 

As in that problem, at the transition point there are 
new critical exponents. More precisely, the Flory expo- 
nent V is the same as for non- grafted animals, but the 
entropic exponent 6 is changed [78]. Since this exponent 
could not be measured by any previous simulation algo- 
rithm and since there exits no field-theoretic predictions 



The critical exponents for animals in d = 2 dimensions 
can be calculated exactly, as for many other critical phe- 
nomena in d = 2 dimensions. But while this is due to 
conformal invariance in these other cases, 2-d animals are 
not conformally invariant [86|. 

For conformally invariant problems of cluster growth, 
the entropic critical exponents of clusters grafted to the 
tips of wedges and cones (wedges with identified edges) 
can be calculated exactly for any wedge angle, by map- 
ping the wedge onto the half plane. Due to lack of confor- 
mal invariance, this is no longer true for 2-d lattice ani- 
mals. In [83], the exponents 9{a) were measured carefully 
not only for wedges and cones with angles up to a = 27r. 
By grafting them to branch points of Riemann sheets, 
angles up to IOtt were studied. Results are shown in Fig. 
1231 The simulations were made with the hope that some- 
one might produce a fit to these data that could suggest 
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FIG. 23. Entropic critical exponents 9{a) for 2-d lattice ani- 
mals grafted to the tips of wedges resp. cones with angles a. 
Adapted from Ref. ]Sl\ . 



an alternative to confornial invariance. So far this hope 
has not materialized, in contrast to what happened 111 
years ago to some obscure black body radiation data [8a| ■ 



D. Collapsing Lattice Animals and Lattice Trees in 

d = 2 

A coil-globule transition similar to that for linear poly- 
mers is also expected to occur for randomly branched 
polymers as the solvent quality becomes worse, but the 
situation is much more complicated. To describe the pos- 
sible collapse transitions for self-interacting lattice an- 
imals, we need two different types of interactions be- 
tween nearest-neighbor monomer-monomer pairs: (cova- 
lent) bonds that are needed for the connectedness of the 
cluster but that can also form loops when present in ex- 
cess, and weak interactions between non-bonded pairs 
("contacts"). Associated to these are two different con- 
trol parameters |89l| . The partition sum is therefore writ- 
ten as follows [90| 



ZN{y,T) 



b,k 



Nbk.y 



-N+1 k 



(61) 



where C^bk is just the number of configurations (up to 
translations and rotations) of connected clusters with 
N sites, b bonds, and k contacts, y and r are fugac- 
ities for monomer-monomer bonds and for non-bonded 
monomer- monomer contacts, respectively. As for un- 
branched polymers, there is no need to introduce a sep- 
arate monomer- solvent interaction, since the number s 
of monomer-solvent contacts is not independent, but is 
given by 



shown in Fig. [24l This model includes the following spe- 
cial cases: 

• Unweighted animals: y — t ^ 1 

• Bond percolation: y — p/{\—pY and r — 1/(1 —p) 
where < p < 1. The critical percolation point is 
at y = 2 and r = 2 as p = pc = 1/2. 

• Collapsing trees: y = where h ^ N — 1. 

• 'Strongly embeddable' animals with t = 0, which 
have no contacts (fc = 0) This model was first stud- 
ied by Derrida and Herrmann by transfer matrix 
methods ^. 

The transition points are determined by the scaling 
laws of the partition sum: 



ZN{y,T^T,)^^i{yY''N 



and the gyration radius 
R 



N 



N" 



(63) 



(64) 



where /i(y) should depend continuously on y, but 6 
should take discrete values depending on the respective 
universality, v is the Flory exponent. A phase diagram 
for interacting animals in d = 2 is shown in Fig. 1251 
Lattice animals are in the extended phase below the full 
line, but they are in the collapsed phase above the full 
line. The bond percolation is described by the dashed 
curve. The percolation critical point at y = 2 and r = 2 
divides the transition line into two different universality 
classes. On the left-hand side, the collapse transitions 
are dominated by non-bonded contacts. In this region 
PERM simulations are very easy and yield very precise 
values for the transition curve (which seems to be ex- 
actly horizontal) and for the critical exponents. These 
results have been fully confirmed by field theoretic meth- 
ods [92]. For the Derrida-Herrmann model at the far 
right end of the transition curve, PERM simulations are 
least efficient, and they could not improve on the results 
of yV\. It is not entirely clear whether there exists a 
further (multi-)critical point between this end point and 
the percolation point. Such a point, together with an 
additional phase separation line emanating from it, was 
suggested by earlier exact enumeration studies (cited in 
(90]). PERM simulations also weakly suggested such an 
addition phase separation line, indicated by the short 
dashed-dotted line in Fig. [551 but these simulations were 
not easy and interpreting their results was not unam- 
biguous. Indeed, a completely different scenario for the 
behavior along the transition line between the percola- 
tion and Derrida-Herrmann points is suggested in [94] . 



NN = 2b + 2k + s , 



(62) 



where M is the lattice coordination number (TV = 2d 
on a simple hypercubic lattice in d = 2 dimensions). 
A schematic drawing of an interacting lattice animal is 



X. PROTEIN FOLDING 

In this section we shall only describe applications of a 
variant of PERM [9J, |93| to simple lattice models, where 
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FIG. 24. A schematic drawing of an interacting lattice animals which contain a cluster with b — 11 bonds and A'^ = 12 sites, 
and k — 2 non-bonded monomer- monomer contacts, and s — 22. It leads to 4A'^ = 2b + 2k + s. 
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FIG. 25. Phase diagram for interacting animals in d = 2. The full curve separates an extended phase (below) from a collapsed 
phase (above). At y — the clusters are trees (minimal number of bonds), while at r = they have no contacts but only 
bonds. The dashed line corresponds to bond percolation, with the critical point being at y = r = 2. The short dashed-dotted 
line is a rough estimate for a possible transition between a contact-rich and a bond-rich collapsed phase. The critical exponents 
u and 6 are also shown for different universality classes. Adapted from Ref. [9Cl |. 



it seems one of the most efRcient algorithms for finding 
low energy states. PERM was also applied to continuum 
models [95] and was there more efficient than previous 
Monte Carlo algorithms |9S] , but has been rendered later 
obsolete in this application [97| . 



A. New Version of PERM (nPERM) 

The main improvement of nPERM is that we no longer 
make identical clones as implemented in old PERM, in 
order to avoid the loss of diversity which limited the suc- 
cess of old PERM. 

When we have a configuration of polymer chains with 
n — 1 monomers, we first estimate a predicted weight 
ly^rod jTqj, ^Yie next step (the n"^ step), and we count 



the number fcfree of free sites where the n monomer can 
be placed. If Wp™'^ > W+ and k^cc > 1, we make k 
{2 < k < fcfroo) clones with the request that k different 
sites are chosen for the n"^ step. Therefore, k configu- 
rations with n monomers are forced to be different. If 
M^P''°'^ < W,7j £^ random number r is chosen uniformly 
in [0,1]. If r < 1/2, the chain is discarded, otherwise 
it is kept and its weight is doubled. We tried several 
strategies for selecting k which all gave similar results. 
Typically, we used k = min {khee, [Wi''"^ /W+]} . 



It is still important to keep the right weight of each con- 
figuration with n monomers. When selecting a /c-tuple 
A — {ai, . . . ,ak} of mutually different continuations Uj 
with probability pA, the corresponding weights Wn,ai, 
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• , Wn,ai, are 



^"."j = ^-TITZT ' J = 1.2,...,fc. (65) 



Here, the importance 



exp{~l3En,aj) 



(66) 



of choice aj is the Bohzmann-Gibbs factor associated 
with the energy En, a of the n*^ placed monomer in the 
potential created by all previous monomers. The other 
terms arise from correcting bias and normalization. 

Two strategies for the choice of k continuations among 
fcfree are described as follows: 

(i) New PERM with simple sampling (nPERMss): 
k different free sites are chosen randomly and uni- 
formly. The predicted weight is given by, 



Wr" = Wn-lh,, 



(67) 



and the corresponding weight for each continuation 



a^ IS 



Wn,a, = Wn-iqajh^cc/k 



(68) 



since there are ( 'jj,'"') different ways to select a k- 
tuple with equal probability, the probability pA is 
therefore 



PA 



Kfr, 



(69) 



Here the tuples related by permutations are con- 
sidered as identical. 

(ii) New PERM with importance sampling (nPERMis): 
k different free sites are chosen according to the 
modified Boltzmann weight g^ . defined by 






+ l/2)exp(-/3K.a,) 



(70) 



where k\^^^ is the number of free neighbors when 
the v}^ monomer is placed at a^, and En^a- is its 
energy gain. The idea of replacing qa- by qa^ is 
that we anticipate continuations with fewer free 
neighbors which will contribute less on the long 
run than continuations with more free neighbors. 
This is similar to " Markovian anticipation" within 
the framework of old PERM, described in section 
Sect. IHIBl where the bias for placing a monomer 
at the next step different from the short-sighted op- 
timal importance sampling was found to be prefer- 
able. The predicted weight is now 



wr 



"-free 

Wn-l Y^ qa, 



(71) 



Using the requirement that the variance of the 
weights Wn is minimal, the proper choice of the 
probability pA to select a tuple A — {ai, . . . , a^} is 
found to be 



PA 






j ^ l,2,...,kh, 



(72) 



If qa had not been replaced by g^ , the variance 
of Wn for fixed Wn-i would be zero. For fc = 1, it 
corresponds to the standard importance sampling, 
i.e. PA = Paj = qaJY!i={ qa,- The weight at the 
n}^ step is thus Wn,aj = Wn-iqaJPa^- For fc > 1, 
Wn,aj is given by ([65]). 

A noteworthy feature of both nPERMss and nPERMis 
is that they cross over to complete enumeration when W^ 
and W~ tend to zero. In this limit, all possible branches 
are followed and none is pruned as long as its weight 
is not strictly zero. In contrast to this, with the use of 
the original PERM as explained in Sect.lHl exponentially 
many copies of the same configuration would be made. 
It suggests that one can be more lenient in choosing W^ 
and Wn when applying nPERM. 



B. HP Model 

For testing the efficiency of the new PERM, we ap- 
plied it to the HP model [Qa] since this model is well 
simulated for bench-marking. In this model, a protein is 
simplified by replacing amino acids by only two types of 
monomers, H (hydrophobic) and P (polar) monomers. 
Therefore a protein (a polymer) of length n is modeled 
as a self-avoiding chain of n steps on a regular (square 
or simple cubic) lattice with repulsive or attractive inter- 
actions between neighboring non-bonded monomers such 
that tHH = — 1, ^HP = ^pp = 0. The partition sum is 



E 



Zn = 2^ q" 

walks 



(73) 



where q = exp(— /3e/f//) and m is the total number of 
non-bonded H — H pairs. 

In our simulations, we chose the two thresholds W„— — 
and Wn < oo, i.e. we neither pruned nor branched, for 
the first configuration hitting length n. For the follow- 
ing configurations we used Wn — CZn/Zo{cn/co)'^ and 
Wn = 0.2 Wn- Here, c„ is the total number of config- 
urations of length n already created during the run, Z„ 
is the partition sum estimated from these configurations, 
and C is some positive number < 1. The idea of incorpo- 
rating the term (c„/co)^ is that we can reduce the upper 
threshold Wn in order to make more cloning in possible 
branches as a lower energy state is hit but only few con- 
figurations of length n have been obtained. The following 
results were all obtained with C = 1, though substantial 
speed-ups (up to a factor 2) could be obtained by choos- 
ing C much smaller, typically as small as 10~^^ to 10~^*. 
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TABLE I. Perfor man ces for the 3-d binary (HP)- sequences 
of 48-mers from 100], presented by the CPU time (minutes) 
per independent ground state hit. The ground state e nerg y 
is denoted by Emin- Results obtained by using PERM [l02l ]. 
nPERMss and nPERMis [H, HI] are carried out on a 167 
MHz Sun Ultra I workstation. Results quoted from Ref. [92| 
obtained by CG are multiplied by 10 for the comparison. 



sequence nr. 


— £'min 


PERM nPERMss nPERMis 


CG 


1 


32 


6.9 


0.66 


0.63 


0.94 


2 


34 


40.5 


4.79 


3.89 


3.50 


3 


34 


100.2 


3.94 


1.99 


6.20 


4 


33 


284.0 


19.51 


13.45 


2.90 


5 


32 


74.7 


6.88 


5.08 


1.20 


6 


32 


59.2 


9.48 


6.60 


46.00 


7 


32 


144.7 


7.65 


5.37 


6.40 


8 


31 


26.6 


2.92 


2.17 


3.80 


9 


34 


1420.0 


378.64 


41.41 


- 


10 


33 


18.3 


0.89 


0.47 


0.11 



The latter is easily understandable: with such small C, 
the algorithm performs essentially exact enumeration for 
short chains, giving thus maximal diversity, and becomes 
stochastic only later when following all possible configu- 
rations would become unfeasible. 

For presenting the efficiency of nPERMss and nPER- 
Mis, we applied them to find the ground state of the HP 
model with blind search. Special comparison is made 
with the core-directed growth method (CG) of Beutler and 
Dill [99|]. This is the only method we found to be still 
competitive with nPERM but it works only for the HP 
model and relies heavily on heuristics. Two examples are 
shown here. 

(a) Ten sequences of 48-mers in d = 3 from Ref. (lO0J | 
are tested. In Table HI we list the required CPU time 
(measured in minutes) per independent ground state hit 
on a 167 MHz Sun ULTRA I workstation. As with the 
original PERM [l4|, we could also reach lowest energy 
states by using nPERM, but the required CPU time is 
within one order of magnitude shorter than that needed 
for PERM. For all ten sequences we use the same tem- 
perature, exp(l/T) = 18, although we could have op- 
timized CPU times by using different temperatures for 
each chain. Results obtained in Ref. [93, l9J] are carried 
out on a SPARC 1 machine which is slower by a factor 
« 10 than the Sun workstation. Therefore, in Table |T] 
we multiplied their results by 10 for comparison. We see 
that nPERM gave comparable speeds as CG. But, one 
has to note that the lowest energy of the sequence No. 9 
was not hit by CG [g^. 

(b) There e xists one HP sequence of 64-mers intro- 
duced in Ref. J10l| . for which it is particularly difficult 
to find its ground state energy £',„in = —42 by any chain 
growth algorithm. One of the ground state configuration 
is shown in Fig. I^STa). Its degeneracies of -Emin = —42 



differ in the detailed folding of the tails in the interior. 
As one uses a chain growth algorithm, it seems very un- 
natural that the chain has to grow first along an arc 
(Fig. BSffb). only until much later that the structure 
of the chain will be stabilized. It shows the difficulty 
of folding this HP sequence into its ground state. With 
nPERM, the ground state was reached with blind search. 
The average CPU time per ground state hit was about 
30h on the DEC21264, which seems to be roughly com- 
parable to the CPU time needed in Refs. (l0lll04l|. but 
slower than Ref. [90| where CG was used. In a previ- 
ous application of (old) PERM [l02| . the configuration 
wirh £^„iin = —42 was found only by means of some spe- 
cial tricks (non-blind search) together with the original 
PERM. 



XI. DNA MELTING 

At physiological temperatures, DNA forms the famous 
double helix. When the temperature is elevated, a point 
Tm is reached where the covalent bonds along the back- 
bone are still strong enough to keep the two strands in- 
tact, but the hydrogen bonds between the strands no 
longer can keep them together. The ensuing separation 
at temperatures > Tm is known as DNA denaturation or 
DNA melting. 

Since the strengths of the hydrogen bonds between 
pairs A-T and C-G are different, also the melting tem- 
peratures Tjn are different for homogeneous DNA, with 
Tm{A -T) <Tm{C -G). For natural DNA, the effec- 
tive melting temperature depends on the A/C compo- 
sition, and precise measurements of melting curves for 
short pieces of DNA can give detailed information about 
the base composition. This has been used for a long time 
as one of the easiest and fastest methods to obtain genetic 
information, and modern developments have made high 
resolution DNA melting one of the most simple, cheap 
and fast techniques for geno typing, sequence matching, 
and mutation scanning jl06| |. 

The sharpness of the transition in case of long homo- 
geneous DNA has suggested since long ago that DNA 
melting is a fir st or der phase transition [105]. But the 
earliest models |l07l | by Poland and Scheraga could only 
give rise to a second order transition, which was seen as 
a severe problem. These models of course lacked many 
aspects of the real DNA melting problem, such as the 
helical structure of DNA. This was done in view of the 
universality of second order phase transitions, and later 
models that did include the helix structure indeed did 
not do better. 

On the other hand, it was already speculated early on 

that the excluded volume effect - that was neglected in 

lOTj - could be responsible for the change into a first 

order transition. The first model that treated the ex- 



cluded volume effect correctly was published by Coluzzi 
et al. jlOSJ I and simulated by means of PERM. The model 
treated each DNA strand as a SAW on the simple cubic 



25 



(a) 



(b) 




O H 
• P 




FIG. 26. (a) One of the ground state configurations for a A'' = 64 cliain in 2D from [101[. Otlier states with the same energy 
differ in the detailed folding of the tails in the interior, but have identical outer shapes, (b) When about 3/4 of the chain is 
grown, one has to pass through a very unstable configuration which is stabilized only later, when the hydrophobic core is filled. 
Adapted from Ref. [H, HJ 



lattice. But the two strands were mutually self avoiding 
only to the extent that bases that were not supposed to 
be bound by hydrogen bonds were not allowed to occupy 
the same lattice site. Base pairs that were bound in the 
native (non-molten) configuration were not only allowed 
to occupy the same site, but would also gain an energy e if 
they did, mimicking thereby the binding between the two 
strands. In addition, variants were studied where either 
the excluded volume effect within each strand and/or be- 
tween the strands was neglected. 

The results were as expected: While all variants 
that did not incorporate the full excluded volume effect 
showed second order transitions, the version with full ex- 
cluded volume interactions showed a first order transi- 
tion. Later studies, both by simulations (using PERM 
and other method s) an d by analytic arguments confirmed 
these results (see |l09l | for references). 



XII. SUMMARY 

In this review we have concentrated on applications 
of PERM to problems in polymer physics. But PERM 
can also be applied to other problems where it is impor- 
tant not only to find rare events, but also to estimate 
the probabilities with which they occur. This includes 
various reaction-diffusion problems such as the long time 
tails in the Donsker-Varadhan problem (see Sect. 3) and 
in the annihilation reaction A + A ^ [llOJ ] , but also to 
more exotic problems like that of multiple spanning clus- 
ters in percolation [18!|. These are all problems where 
theory makes clear predictions that were very hard to 
verify numerically with other algorithms. But there are 
also some other applications of PERM to polymer prob- 
lems that we have not discussed here, such as polymers 
grafted to porous [llll | and non-porou s [5211 membranes. 



adsorption of copolymers to surfaces 
rections for SAWs on Manhattan lattices 
ISAWs with orientation dependent interactions 
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PERM belongs to a class of Monte Carlo algorithms 



called some times "sequential algorithms with resam- 
pling" [ll5|. In contrast to most other algorithms in 
this class, it is implemented depth-first which leads to 
very compact codes and minimal memory requirements. 
In principle, it can be applied to any problem where in- 
stances are built sequentially by repeating small steps. 
Its main ideas are that these steps can be biased in order 
to shear the evolution towards the wanted (in general 
rare but highly weighted) configurations. If this is not 
deemed successful, the further evolution can be pruned, 
while very successful trials can be cloned, with each clone 
evolving further independently. Notice that pruning and 
cloning are done on partially constructed configurations, 
with the hope that configurations that are successful at 
an early stage will also continue to be successful later. 
When this is true, efficiency can be spectacular (such as 
for <d polymers. Sect. 4). But when it is not true, the 
method simply fails. Examples of the latter were also 
discussed in Sects. 3, 4 and 11. 

In most applications, the criterion for success is simply 
the weight of the configuration, based on a combination 
of Boltzmann, entropic, and bias compensating factors. 
But in some cases - illustrated in Sect. 9 for lattice an- 
imals - the weight itself would be a very poor "fitness" 
indicator. For lattice animals, a much better fitness func- 
tion was found empirically. 

In addition to the versions of PERM that we hav e dis- 
cussed in this re view , there exist also "flat" |ll6l | and 
"multicanonical" |ll7l | versions of it. Their main advan- 
tage is that data over a wide range of energies can be 
obtained in one single run, while ordinary PERM would 
need several runs, each covering the energy range that 
dominates at one particular temperature. This is cer- 
tainly an attractive feature, but it is not as important 
as in Markov Chain Monte Carlo algorithms. While it 
is there highly non-trivial to combine results obtained at 
different temperatures [32| , this is much easier for PERM 
where the algorithm provides very precise estimates of 
the partition sum. 
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